Ajuste
INFORME SOBRE AJUSTE NO LINEAL Autores: Carlos Mohor, Pablo Rios, Alfredo Valenzuela Curso: Métodos Numéricos Profesor: Stefan Berres Universidad Católica de Temuco Temuco, 27 de Noviembre del 2012 'Resumen Ejecutivo' El presente informe expone el tema de análisis de modelos no-lineales. El objetivo es que, dado un modelo no-lineal, con variables desconocidas, podamos encontrar los parámetros que ajusten el modelo de mejor manera a los datos observados. Una vez que podamos entender bien la lógica, la idea es poder formular el algoritmo para implementarlo en código y analizar los resultados. Para esto hemos estructurado el informe de la siguiente manera: I. Introducción II. Explicación paso a paso de cómo buscar los parámetros, linealizar el modelo y obtener las ecuaciones en base a un modelo de función logística trabajado en clases. III. Formulación de un algoritmo general en MATLAB IV. Implementación paso a paso de un algoritmo para un modelo no-lineal distinto (variante). V. Implementación paso a paso del algoritmo pero en un modelo nuevo (no visto en clases), previa introducción. 'I. Introducción' ¿Qué es un modelo matemático? • Conjunto de ecuaciones que relacionan las variables de interés del proceso y representan adecuadamente su comportamiento • Siempre son aproximaciones de la realidad • Distintos modelos para distintos objetivos y tipos de procesos • Compromiso entre facilidad de uso y exactitud "NO LINEALIDAD" ' Los sistemas no lineales representan sistemas cuyo comportamiento no es expresable como la suma de los comportamientos de sus descriptores. Más formalmente, un sistema físico, matemático o de otro tipo es no lineal cuando las ecuaciones de movimiento, evolución o comportamiento que regulan su comportamiento son no lineales. En particular, el comportamiento de sistemas no lineales no está sujeto al principio de superposición, como lo es un sistema lineal. La linealidad de un sistema permite a los investigadores hacer ciertas suposiciones matemáticas y aproximaciones, permitiendo un cálculo más sencillo de los resultados. Ya que los sistemas no lineales no son iguales a la suma de sus partes, usualmente son difíciles (o imposibles) de modelar, y sus comportamientos con respecto a una variable dada (por ejemplo, el tiempo) es extremadamente difícil de predecir. Algunos sistemas no lineales tienen soluciones exactas o integrables, mientras que otros tienen comportamiento caótico, por lo tanto no se pueden reducir a una forma simple ni se pueden resolver. Un ejemplo de comportamiento caótico son las olas gigantes. Aunque algunos sistemas no lineales y ecuaciones de interés general han sido extensamente estudiados, la vasta mayoría son pobremente comprendidos. 'II. Explicación paso a paso el desarrollo del modelo de función logística OBJETIVO: Entender la lógica del desarrollo dado el modelo: f(x)=\frac{a}{b+c\exp(-dx)} y los datos \begin{array} \hline X_i & 1 & 2 & 3 & 4 & 5 & 6\\ \hline f_i & 1 & 4 & 15 & 30 & 50 & 55\\ \hline \end{array} formular el algoritmo que minimiza S= \sum_{i=1}^n f_i-f(x_i)^2 buscar los parametros a, b, c, d que minimiza S Criterio de extreme(minimo) \frac{ds}{da} = 0 \frac{ds}{db} = 0 \frac{ds}{dc} = 0 \frac{ds}{dd} = 0 Tenemos 4 parametros (como variables desconocidas) y 6 ecuaciones(que corresponden a (x_i,f_i), i=1, ... , 6) ) Linealizar con respecto a \begin{array}{llll}\Delta a, \Delta b, \Delta c, \Delta d \end{array} Buscar y calcular los parametros a^*, b^*, c^*, d^* , tal que S(a^*, b^*, c^*, d^*)\min S(a, b, c, d) Parametros a, b, c, d para que la curva (x,f(x)) se ajuste de mejor manera a los datos. Los parametros deben coincidir con la escala, se ajustan a la cuerva(podrian tener la misma forma, pero en otra escala). Paso1: "Transformar eñ problema de minimizacion a una tarea de calcular ceros" Queremos minimizar S(a, b, c, d) \begin{array}{llll} \left . \begin{matrix} F_1(a,b,c,d):=\frac{\partial S(a,b,c,d)}{\partial a}=0\\ F_2(a,b,c,d):=\frac{\partial S(a,b,c,d)}{\partial b}=0\\ F_3(a,b,c,d):=\frac{\partial S(a,b,c,d)}{\partial c}=0\\ F_4(a,b,c,d):=\frac{\partial S(a,b,c,d)}{\partial d}=0\\ \end{matrix} \right \}=F=\Delta S(a,b,c,d)\end{array} Definir: Todos los componentes de la funcion F dependen de los parametros(dependencia). F = \left ( \begin{matrix} F_1 \\ F_2 \\ F_3 \\ F_4 \end{matrix} \right ) , F = F(a,b,c,d) \begin{array}{llll}F =\Delta S; \end{array} (F corresponde a la gradiente de la funcion S). Para encontrar el minimo de S(a,b,c,d) necesitamos encontrar el cero de F: F(a,b,c,d) =\left ( \begin{matrix} 0 \\ 0 \\ 0 \\ 0 \end{matrix} \right ) Calculando los componentes de F: F_1(a,b,c,d) =\frac{\partial S(a,b,c,d)}{\partial a} = \frac{\partial }{\partial a} \sum_{i=1}^n (f(x_i)-f_i)^2) = \sum_{i=1}^n \frac{\partial }{\partial a}(f(x_i)-f_i)^2) = \sum_{i=1}^n 2(f(x_i)-f_i) \frac{\partial }{\partial a} f(x_i) ... lo mismo para los otros F_2, F_3, F_4 ... Por lo tanto; F_1(a,b,c,d) = \left [ \sum_{i=1}^n 2(f(x_i)-f_i) \frac{\partial }{\partial a} f(x_i) \right ] F_2(a,b,c,d) = \left [ \sum_{i=1}^n 2(f(x_i)-f_i) \frac{\partial }{\partial a} f(x_i) \right ] F_3(a,b,c,d) = \left [ \sum_{i=1}^n 2(f(x_i)-f_i) \frac{\partial }{\partial a} f(x_i) \right ] F_4(a,b,c,d) = \left [ \sum_{i=1}^n 2(f(x_i)-f_i) \frac{\partial }{\partial a} f(x_i) \right ] Paso2: "Linealizar el modelo no lineal" Dificultad: -Entender y aplicar rotacion vectorial. -Entender la secuancia de varios pasos de transformacion del problema Linealizar el modelo no lineal es decir, linealizar: F(a,b,c,d) Para aproximar F(P)=0 , con P=(a,b,c,d) , Introduciremos P_i=(a_i,b_i,c_i,d_i) , i=0,1,2,... Formular el metodo de newton, con la especificacion de la funcion: F(a,b,c,d) F(P_{i+1})=F(p_i)+\left ( \begin{matrix} \frac{\partial F_1}{\partial a} & \frac{\partial F_1}{\partial b} & \frac{\partial F_1}{\partial c} & \frac{\partial F_1}{\partial d}\\ \frac{\partial F_2}{\partial a} & \frac{\partial F_2}{\partial b} & \frac{\partial F_2}{\partial c} & \frac{\partial F_2}{\partial d}\\ \frac{\partial F_3}{\partial a} & \frac{\partial F_3}{\partial b} & \frac{\partial F_3}{\partial c} & \frac{\partial F_3}{\partial d}\\ \frac{\partial F_4}{\partial a} & \frac{\partial F_4}{\partial b} & \frac{\partial F_4}{\partial c} & \frac{\partial F_4}{\partial d} \end{matrix} \right ) \Delta P = 0 Iteracion con respecto a i \to i+1 (1) g(P_i)\Delta P = -F(p_i) \to con esto podremos implementar el algoritmo (2) P_{i+1}=P_i+\Delta P donde \Delta P=P_{i+1}-P_i=\left ( \begin{matrix} a_{i+1}- a_i \\ b_{i+1}- b_i \\ c_{i+1}- c_i \\ d_{i+1}- d_i \end{matrix} \right )=\left ( \begin{matrix} \Delta a \\ \Delta b \\ \Delta c \\ \Delta d \end{matrix} \right ) -Especificar la matriz g(P_i) , en termino de solamente a, b, c, d o sea: -Especificar g(P_i) en terminos de f ,\frac{\partial f}{\partial a},\frac{\partial f}{\partial b},\frac{\partial f}{\partial c},\frac{\partial f}{\partial d} -Calcular derivadas parciales de f . g(P_i) , en terminos de f,a ,b ,c ,d: Tenemos; g(p_i)=\left ( \begin{matrix} g_{11}(p_i) & g_{12}(p_i) & g_{13}(p_i) & g_{14}(p_i)\\ g_{21}(p_i) & g_{22}(p_i) & g_{23}(p_i) & g_{24}(p_i)\\ g_{31}(p_i) & g_{32}(p_i) & g_{33}(p_i) & g_{34}(p_i)\\ g_{41}(p_i) & g_{42}(p_i) & g_{43}(p_i) & g_{44}(p_i) \end{matrix} \right ) \begin{matrix} g_{11}(p_i)\frac{\partial F_1}{\partial a} & g_{12}(p_i)\frac{\partial F_1}{\partial b} & g_{13}(p_i)\frac{\partial F_1}{\partial c} & g_{14}(p_i)\frac{\partial F_1}{\partial c}\\ g_{21}(p_i)\frac{\partial F_1}{\partial a} & g_{22}(p_i)\frac{\partial F_1}{\partial b} & g_{23}(p_i)\frac{\partial F_1}{\partial c} & g_{24}(p_i)\frac{\partial F_1}{\partial c}\\ g_{31}(p_i)\frac{\partial F_1}{\partial a} & g_{32}(p_i)\frac{\partial F_1}{\partial b} & g_{33}(p_i)\frac{\partial F_1}{\partial c} & g_{34}(p_i)\frac{\partial F_1}{\partial c}\\ g_{41}(p_i)\frac{\partial F_1}{\partial a} & g_{42}(p_i)\frac{\partial F_1}{\partial b} & g_{43}(p_i)\frac{\partial F_1}{\partial c} & g_{44}(p_i)\frac{\partial F_1}{\partial c} \end{matrix} \to calculamos: (1) Para g_{11}(P_i) ; g_{11}(P_i) = \frac{\partial F_1}{\partial a}; \to se aplica la definicion de: F_1=f_1(a,b,c,d):= \frac{\partial S}{\partial a}(a,b,c,d) \therefore g_{11}(P_i)= \frac{\partial F_1}{\partial a} =\frac{\partial }{\partial a}\left [ \sum_{i=1}^n 2(f(x_i)-f_i) \frac{\partial }{\partial a} f(x_i) \right ] = \sum_{i=1}^n \frac{\partial }{\partial a}\left \frac{\partial }{\partial a} f(x_i) \right =\sum_{i=1}^n \left }{\partial a}\left (\frac{\partial }{\partial a}f(x_i)\right )+\frac{\partial }{\partial a}f(x_i)\frac{\partial }{\partial a} 2(f(x_i)-f_i)\right =\left ^2}{\partial a^2}f(x_i) + 2\left (\frac{\partial }{\partial a}f(x_i)\right )^2\right g_{11}(P_i) =\sum_{i=1}^n 2\left ^2}{\partial a^2}f(x_i) + \left (\frac{\partial }{\partial a}f(x_i)\right )^2\right (2) Para g_{12}(P_i) ; \frac{\partial F_1}{\partial b}=2f(x_i)\frac{\partial }{\partial b}\left (\frac{\partial }{\partial a}f(x_i)\right )+\frac{\partial }{\partial a} f(x_i)\frac{\partial }{\partial b}2f(x_i) (3) Para g_{13}(P_i) ; \frac{\partial F_1}{\partial c}=2f(x_i)\frac{\partial }{\partial c}\left (\frac{\partial }{\partial a}f(x_i)\right )+\frac{\partial }{\partial a} f(x_i)\frac{\partial }{\partial c}2f(x_i) (4) Para g_{14}(P_i) ; \frac{\partial F_1}{\partial d}=2f(x_i)\frac{\partial }{\partial d}\left (\frac{\partial }{\partial a}f(x_i)\right )+\frac{\partial }{\partial a} f(x_i)\frac{\partial }{\partial d}2f(x_i) (5) Para g_{21}(P_i) ; \frac{\partial F_2}{\partial a}=2f(x_i)\frac{\partial }{\partial a}\left (\frac{\partial }{\partial b}f(x_i)\right )+\frac{\partial }{\partial b} f(x_i)\frac{\partial }{\partial a}2f(x_i) (6) Para g_{22}(P_i) ; \frac{\partial F_2}{\partial b}=2\left ^2}{\partial b^2}f(x_i) + \left (\frac{\partial }{\partial b}f(x_i)\right )^2\right (7) Para g_{23}(P_i) ; \frac{\partial F_2}{\partial c}=2f(x_i)\frac{\partial }{\partial c}\left (\frac{\partial }{\partial b}f(x_i)\right )+\frac{\partial }{\partial b} f(x_i)\frac{\partial }{\partial c}2f(x_i) (8) Para g_{24}(P_i) ; \frac{\partial F_2}{\partial d}=2f(x_i)\frac{\partial }{\partial d}\left (\frac{\partial }{\partial b}f(x_i)\right )+\frac{\partial }{\partial b} f(x_i)\frac{\partial }{\partial d}2f(x_i) (9) Para g_{31}(P_i) ; \frac{\partial F_3}{\partial a}=2f(x_i)\frac{\partial }{\partial a}\left (\frac{\partial }{\partial c}f(x_i)\right )+\frac{\partial }{\partial c} f(x_i)\frac{\partial }{\partial a}2f(x_i) (10) Para g_{32}(P_i) ; \frac{\partial F_3}{\partial b}=2f(x_i)\frac{\partial }{\partial b}\left (\frac{\partial }{\partial c}f(x_i)\right )+\frac{\partial }{\partial c} f(x_i)\frac{\partial }{\partial b}2f(x_i) (11) Para g_{33}(P_i) ; \frac{\partial F_3}{\partial c}=2\left ^2}{\partial c^2}f(x_i) + \left (\frac{\partial }{\partial c}f(x_i)\right )^2\right (12) Para g_{34}(P_i) ; \frac{\partial F_3}{\partial d}=2f(x_i)\frac{\partial }{\partial d}\left (\frac{\partial }{\partial c}f(x_i)\right )+\frac{\partial }{\partial c} f(x_i)\frac{\partial }{\partial d}2f(x_i) (13) Para g_{41}(P_i) ; \frac{\partial F_4}{\partial a}=2f(x_i)\frac{\partial }{\partial a}\left (\frac{\partial }{\partial d}f(x_i)\right )+\frac{\partial }{\partial d} f(x_i)\frac{\partial }{\partial a}2f(x_i) (14) Para g_{42}(P_i) ; \frac{\partial F_4}{\partial b}=2f(x_i)\frac{\partial }{\partial b}\left (\frac{\partial }{\partial d}f(x_i)\right )+\frac{\partial }{\partial d} f(x_i)\frac{\partial }{\partial b}2f(x_i) (15) Para g_{43}(P_i) ; \frac{\partial F_4}{\partial c}=2f(x_i)\frac{\partial }{\partial c}\left (\frac{\partial }{\partial d}f(x_i)\right )+\frac{\partial }{\partial d} f(x_i)\frac{\partial }{\partial c}2f(x_i) (16) Para g_{44}(P_i) ; \frac{\partial F_4}{\partial c}=2\left ^2}{\partial d^2}f(x_i) + \left (\frac{\partial }{\partial d}f(x_i)\right )^2\right 'III. Formulación de un algoritmo general en MATLAB' Seudocódigo para una implementación en MATLAB: Basado en el modelo f(a,b;x) Sumatoria “S” es la misma que hemos utilizado hasta ahora, con los datos Xi,Fi %Declaracion de variables %generacion de datos (x(i),f(i)) a=1, b=2 for i=1 hasta n %generamos los datos x(i) = i/n; f(i) = a*b*x(i) end for a=0, b=1 %escogemos los parametros iniciales de iteración i=0; f=1;1; épsilon = 1e-7; xx=[] %Iteracion: While norm(F)<épsilon && i %calculamos F1 (la derivada de S con respect a A) %f(a,b;x)=abx F1=0; For j=1:n , F1 = F1+2*(a*b*x(i)-f(i))*b*x(i) End for J = J12, J21, J22 F = F1,F2 Dx=-J/F %Resolver sistema de ecuaciones X=X+Dx %a=x(1);b=X(2); XX=XX,X End while %visualizacion %Plot() 'IV. Implementación paso a paso de un algoritmo para un modelo no-lineal distinto (variante). Implementado en código' ____________________________________________________________________________________ I. PRIMERA IMPLEMENTACION: La primera implementacion esta hecha en base a la aplicacion de un modelo variante expuesto en la clase del martes 06 de noviembre del 2012. Este modelo busca continuar experimentando con el trabajo de minimizacion. El modelo corresponde a un f(t) = A*exp(B*t) con parametros A, B y los datos (ti,fi), i = 1,2,3, ... n. Se busca formular el algoritmo que minimiza la sumatoria trabajada anteriormente. ____________________________________________________________________________________ ______________________________________________________________________________ SOBRE EL CODIGO: Este codigo esta explicado en detallo de la misma forma que explicamos al principio los pasos a seguir para desarrollar un ajuste de datos a un modelo-no lineal. En el caso anterior fue con el modelo de funcion logistica. Para esta primera implementacion utilizaremos otro modelo. El codigo esta implementado en python con una libreria matematica que nos permite resolver los problemas y que trabaja en terminos muy similares a MATLAB. SYMPY: SymPy es una biblioteca escrita en Python cuyo objetivo es reunir todas las caracteristicas de un sistema de algebra computacional (CAS), ser facilmente extensible y mantener el codigo todo lo simple que sea posible. SymPy no requiere ninguna biblioteca externa, salvo para soporte grafico. Caracteristicas: En su funcionalidad podemos distinguir entre: Capacidades basicas, que incluyen:manejo de enteros de precision arbitraria y de numeros racionales, simplificacion basica, expansion, sustitucion basica,manejo de funciones sobre el cuerpo de los complejos, derivacion, expansion en series de Taylor o de Laurent,simbolos no conmutativos. Modulos que incorporan estas tareas: mas funciones (factorial, zeta, legendre, etc),limites,integracion,divisibilidad y factorizacion de polinomios, resolucion de ecuaciones algebraicas, diferenciales y sistemas,operaciones con matrices simbolicas, Algebra de Dirac y de Pauli, Representacion grafica (en 2D y en 3D). O paquetes externos: symbide: GUI en PyGTK --> mas informacion sobre sympy en el siguiente link: http://docs.sympy.org/0.7.2/modules/solvers/solvers.html# ____________________________________________________________________________________ ''--> CÓDIGO!'' #importamos las librerias de sympy... para calculos y variables que #trabajamos como symbolos import sympy from sympy.abc import x, a, b, c import scipy from sympy import summation, oo, symbols, log, exp from sympy import Symbol from sympy import * from sympy.plotting import plot from sympy.plotting import plot3d #___________________DEFINIMOS EL MODELO (FUNCION f) A = Symbol('A') B = Symbol('B') t = Symbol('t') f = Function("f") expresion = A*exp(B*t) #--> MODELO! f = expresion #print f #print f.diff(t) #para derivar #print diff(f,t) #para derivar #f = A**2 #A=2 #probando reemplazos #f = A**2 #probando reemplazos #print "f(2) = " + str(f) #probando reemplazos #g = Function("g") #g = t**2 #print g.subs(t,4).evalf() #para evaluar la funcion f, con respecto a t en un valor cualquiera #____________#evaluando la funcion f(t), con t = 4 #print f.subs(t,4).evalf() #__________________________________________________ #_________________DEFINIMOS LA SUMATORIA S ti = Symbol('ti') fi = Symbol('fi') i = Symbol('i') #para iniciar suma n = Symbol('n') #para termino suma z = Symbol('z') #funcion a evaluar #z = (f(ti)-fi)**2 --> expresion de la sumatoria z = (f.subs(t,ti).evalf() - fi)**2 #expresion de la sumatoria S, reemplazada en la funcion F # con expresion = A*exp(B*t) --> f(ti) #n = 2 #S = summation(z, (i, 1, n)) #sumatoria de S: desde i=1 hasta n #print S ## PASO1: "Transformar el problema de minimizacion a una tarea de calcular ceros" # Queremos minimzar S(A,B) #Calculamos los componentes del vector F: F1 = Symbol('F1') F2 = Symbol('F2') F1 = diff(z,A) F2 = diff(z,B) #print F1 #print F2 #Una vez que tenemos F1 y F2 vamos al paso2 ## PASO2: "Linealizar el modelo no-lineal" #Dificultad: Entender y aplicar notacion vectorial #Entender la secuencia de varios passos de transformacion del problema #Linealizar el modelo no-lineal, es decir, linealizar #F(A,B) #para aproximar vector F(P) = 0, con vector P=(A,B) #Introducimos Pi= (Ai, Bi), con i=0,1,2... #formulamos el metodo de newton con la especificacion de F(A,B) # J11 = Diferencial de F1 con respcto a A # J12 = Diferencial de F1 con respecto a B # J21 = Diferencial de F2 con respecto a A # J22 = Diferencial de F2 con respecto a B J11 = Symbol('J11') J12 = Symbol('J12') J21 = Symbol('J21') J22 = Symbol('J22') J11 = diff(F1,A) #term.B #print J11 J12 = diff(F1,B) #term.AyB #print J12 J21 = diff(F2,A) #term.AyB #print J21 J22 = diff(F2,B) #term.AyB #print J22 ecuaciones_J = J11,J12,J21,J22 #print ecuaciones3 epsilon = 0.0001 F = F2 #resolver sistema de ecuaciones ... usamos from sympy.solvers import solve equations = [ Eq(-J11, f), Eq(-J12, f), Eq(-J21, f), Eq(-J22, f), ] #print solve(Eq(-J11, f)) n = 40 valor_ti = 2 #podemso jugar con los parametros valor_fi = 4 #en distintos valores result = [] J11 = J11.subs(ti,valor_ti).evalf() sol1 = summation(J11, (i, 1, n)) #desde i=1 hasta n aux = solve(-J11) #result.append( aux ) #print aux #print J12 J12 = J12.subs(ti,valor_ti).evalf() J12 = J12.subs(fi,valor_fi).evalf() sol2 = summation(J12, (i, 1, n)) #desde i=1 hasta n print sol2 aux = solve(-J12, A) #cambiar si buscamos A o B result.append( aux0 ) #print result0.subs(B,5).evalf() #evaluando valores J21 = J21.subs(ti,valor_ti).evalf() J21 = J21.subs(fi,valor_fi).evalf() sol2 = summation(J21, (i, 1, n)) #print sol2 aux = solve(-J21, A) result.append( aux0 ) J22 = J22.subs(ti,valor_ti).evalf() J22 = J22.subs(fi,valor_fi).evalf() sol2 = summation(J22, (i, 1, n)) #print sol2 aux = solve(-J22, A) result.append( aux0 ) Imagen del código en el editor: thumb|left|400px|Imagen del código en el editor 'V. Implementación paso a paso del algoritmo pero en un modelo nuevo (no visto en clases).' Modelo Propio Si los parámetros del modelo entran en la ecuación en forma no lineal, entonces tenemos un modelo no lineal ¿Lineal o no lineal? E(Y)=\begin{array}{llll} \beta_0 \ + \beta_1 \ x_1 +...+ \beta_k \ x_k \qquad\qquad E(Y) = \frac{ \beta_0 \ }{ \beta_1 \ + \beta_2 \ x )}^{ \beta_3 \ }} \end{array} E(Y)= \begin{array}{llll} \frac{1}{ \beta_0 \ + \beta_1 \ x + \beta_2 \ x^2 } \end{array} E(Y)= \begin{array}{llll} \beta_0 \ exp(-exp(- \beta_0 \ + \beta_1 \ x)) \qquad E(Y)= \beta_0 \ + \beta_1 \ x+ \beta_2 \ x^2 +...+ \beta_k \ x^k \end{array} Muchas aplicaciones en biología, medicina,química, agricultura. En muchos casos surgen a partir de mecanismos físicos, químicos o biológicos conocidos(usando por ejemplo, ecuaciones diferenciales). Permiten modelar datos reales en forma “natural”con mucha flexibilidad(por ejemplo, con asíntotas,valores positivos de las Y, un único valor máximo). Muchos parámetros tienen interpretación práctica(tasas de degradabilidad, crecimiento máximo,velocidad máxima de absorción, etc.) Problemas con los modelos no lineales Se deben usar métodos iterativos para estimar los parámetros por mínimos cuadrados, máxima verosimilitud, etc. Es común encontrar problemas de convergencia en los métodos iterativos. Distintas parametrizaciones pueden afectar no solo la interpretación sino también las propiedades de los estimadores(por ej., una parametrización puede ser muy interesante desde el punto de vista de su interpretación práctica, pero muy mala para lograr convergencia de estimadores). Familias de funciones no lineales:curvas de crecimiento Exponencial: \begin{array}{llll} \frac{dE(Y)}{dt}= \gamma \ E(y);\quad E(Y) = \alpha \ exp( \gamma \ t) \end{array} Monomolecular: \begin{array}{llll} \frac{dE(Y)}{dt}= \gamma \ ( \beta \ -Y);\quad E(Y) = \beta_0 \ - \alpha \ exp(- \gamma \ t) \end{array} Logístico: \begin{array}{llll} \frac{dE(Y)}{dt}= \gamma \ E(Y)( \beta \ - E(Y));\quad E(Y) = \frac{ \beta \ }{1+exp(- \alpha \ + \gamma \ t)} \end{array} Gompertz: \begin{array}{llll} \frac{dE(Y)}{dt}= \gamma \ E(Y)\beta \ - log E(Y) ;\quad E(Y) = \beta \ exp\alpha \ exp(- \gamma \ t ) \end{array} Richards: \begin{array}{llll} \frac{dE(Y)}{dt}=- \frac{ \gamma \ }{ \kappa \ } E(Y)\frac{E(y)}{ \beta \ })^{ \kappa \ } - 1;\quad E(Y) = \beta \ (1 + \kappa \ exp(- \gamma \ (t- \zeta \ )))^{ \frac{-1}{ \kappa \ }} \end{array} Ejemplos de otras familias de funciones no lineales Modelo lineal con puntos de cambio: \begin{array}{llll} E(Y)= \beta_0 \ + \beta_1 \ x I(x \le \alpha \ ) + ( \beta_1 \ a+ \beta_2 \ (x- \alpha \ ))I(x> \alpha \ ) \end{array} Modelo polinomial recíproco(Holliday): \begin{array}{llll} E(Y)= ( \alpha \ + \beta \ x + \gamma \ x^2)^{-1} \end{array} Dentro de los llamados modelos no lineales, existen 3 tipos; Reciproco, polinomial y logistico. Para la Implementacion numero 2, trabajaremos con un modelo Reciproco. '--> CÓDIGO:' Otro modelo NO-Lineal import sympy from sympy.abc import x, a, b, c import scipy from sympy import summation, oo, symbols, log, exp from sympy import Symbol from sympy import * from sympy.plotting import plot from sympy.plotting import plot3d import pylab as P from matplotlib.transforms import offset_copy #___________________DEFINIMOS EL MODELO (FUNCION f) A = Symbol('A') B = Symbol('B') Y = Symbol('Y') x = Symbol('x') f = Function("f") expresion = (A+B*x+Y*(x**2))**(-1) #--> MODELO polinomial reciproco (Holliday) f = expresion ti = Symbol('ti') fi = Symbol('fi') i = Symbol('i') #para iniciar suma n = Symbol('n') #para termino suma z = Symbol('z') #funcion a evaluar z = (f.subs(Y,ti).evalf() - fi)**2 #expresion de la sumatoria S, reemplazada en la funcion F # con expresion = A*exp(B*t) --> f(ti) F1 = Symbol('F1') F2 = Symbol('F2') F1 = diff(z,A) F2 = diff(z,B) J11 = Symbol('J11') J12 = Symbol('J12') J21 = Symbol('J21') J22 = Symbol('J22') J11 = diff(F1,A) J12 = diff(F1,B) J21 = diff(F2,A) J22 = diff(F2,B) ecuaciones_J = J11,J12,J21,J22 epsilon = 0.0001 F = F2 equations = [ Eq(-J11, f), Eq(-J12, f), Eq(-J21, f), Eq(-J22, f), ] n = 40 valor_ti = 2 #podemso jugar con los parametros valor_fi = 4 #en distintos valores result = [] J11 = J11.subs(ti,valor_ti).evalf() sol1 = summation(J11, (i, 1, n)) #desde i=1 hasta n aux = solve(-J11) J12 = J12.subs(ti,valor_ti).evalf() J12 = J12.subs(fi,valor_fi).evalf() sol2 = summation(J12, (i, 1, n)) #desde i=1 hasta n print sol2 aux = solve(-J12, A) #cambiar si buscamos A o B result.append( aux0 ) J21 = J21.subs(ti,valor_ti).evalf() J21 = J21.subs(fi,valor_fi).evalf() sol2 = summation(J21, (i, 1, n)) aux = solve(-J21, A) result.append( aux0 ) J22 = J22.subs(ti,valor_ti).evalf() J22 = J22.subs(fi,valor_fi).evalf() sol2 = summation(J22, (i, 1, n)) aux = solve(-J22, A) result.append( aux0 ) #EXPERIMENTOS DE VISUALIZACION print str(result) #plot(result0, (B, -5, 5)) #plot3d(J12, (A, -5, 5), (B, -5, 5)) Para la visualizacion vamos a experimentar con la libreria MATPLOTLIB.. esta libreria esta desarrollada para python en base a MATLAB.. los conceptos y sintaxis son muy parecidos... para mas informacion sobre documentacion y tipos de graficos en el siguiente link: http://matplotlib.org/gallery.html import pylab as P from matplotlib.transforms import offset_copy X = P.arange(7) Y = X**2 fig = P.figure(figsize=(5,5)) ax = P.subplot(2,1,1) transOffset = offset_copy(ax.transData, fig=fig, x = (result0.subs(B,7).evalf()), #estamos probando valores y=(result0.subs(B,5).evalf()), units='inches') for x, y in zip(X, Y): P.plot((x,),(y,), 'ro') P.show() Gráfica de resultados: thumb|left|400px